The momentum distribution of the homogeneous electron gas 
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We calculate the off-diagonal density matrix of the homogeneous electron gas at zero tempera- 
ture using unbiased Reptation Monte Carlo for various densities and extrapolate the momentum 
distribution, and the kinetic and potential energies to the thermodynamic limit. Our results on the 
renormalization factor allows us to validate approximate Go Wo calculations concerning quasiparticle 
properties over a broad density region (1 < r s < 10) and show that near the Fermi surface, vertex 
corrections and self-consistency aspects almost cancel each other out. 
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The uniform electron gas (jellium) is one of the most 
fundamental models for understanding electronic proper- 
ties in simple metals and semiconductors. Knowledge of 
its ground state properties, and, in particular, of modifi- 
cations due to electron correlation are at the heart of all 
approximate approaches to the many-electron problem in 
realistic models. Quantum Monte Carlo methods (QMC) 
PQ have provided the most precise estimates of the corre- 
lation energy, electron pair density and structure factor 
of jellium; basic quantities for constructing and parame- 
terizing the exchange-correlation energy used in density 
functional theory (DFT) g]. 

Correlations modify the momentum distribution, n^, 
of electrons, and introduce deviations from the ideal 
Fermi-Dirac step-function. The magnitude of the dis- 
continuity at the Fermi surface {hp), the renormaliza- 
tion factor Z , quantifies the strength of a quasi-particle 
excitation [3] and plays a fundamental role in Fermi liq- 
uid and many-body perturbation theory (GW) for spec- 
tral quantities. Whereas the momentum distribution (as 
well as other spectral information) is inaccessible in cur- 
rent Kohn-Sham DFT formulations, the reduced single- 
particle density matrix - the Fourier transform of rik in 
homogeneous systems - is the basic object in the so-called 
density-matrix functional theory these theories rely 
on knowledge of of jellium. Inelastic x-ray scattering 
measurement of the Compton profile of solid sodium [3] 
have determined nk, but experiments for elements with 
different electronic densities are less conclusive. 

In this paper, we calculate for the electron gas (jel- 
lium) by QMC in the density region 1 < r s < 10. Here, 
r s = (47ma^/3) -3 is the Wigner-Seitz density parame- 
ter, n is the density, and clb — fr 2 /me 2 is the Bohr radius. 
In contrast to previous calculations [6], our calculations 
are based on more precise backflow (BF) wave functions 
[7], and a careful extrapolation to the thermodynamic 
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FIG. 1: The momentum distribution (n*,) of the unpolarizcd 
electron gas for various densities extrapolated to the thermo- 
dynamic limit. The inset shows the extrapolation of nt for 
r a = 5 from a system with N = 54 electrons to the thermo- 
dynamic limit, iV — > oo, leading to a significant reduction of 
the renormalization factor Z. 



limit H ■ Similar to the worm algorithm in finite tem- 
perature path-integral and lattice Monte Carlo [lOj [11] , 
we have extended Reptation Monte Carlo (RMC) to 
include the off-diagonal density matrix in order to ob- 
tain an unbiased estimator of the momentum distribu- 
tion |1 31 114] . From our extrapolation scheme, we derive 
the exact behavior of nt close to the Fermi surface. By 
comparing the renormalization factor, Z, with different 
approximate GW theories, we can judge the importance 
of self-consistency and vertex corrections within these ap- 
proaches. The excellent agreement of our QMC results 
with Go Wo over a broad density region indicate strong 
cancellations of vertex and self-consistency corrections 
close to the Fermi surface. 
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r s 


1 


2 


3.99 


5 


10 


E 


1.173(2) 


0.0039(1) 


-0.1555(1) 


-0.1520(1) 


-0.1071(1) 


T 


2.290(3) 


0.6024(5) 


0.1688(1) 


0.1131(1) 


0.0349(1) 


V 


-1.116(1) 


-0.5985(1) 


-0.3243(1) 


-0.2651(1) 


-0.1421(1) 


g(0) 


0.268(3) 


0.152(2) 


0.057(2) 


0.034(1) 


0.0036(4) 


n 


0.999 


0.998 


0.97 


0.93 


0.88 


n 2 


0.038 


0.066 


0.12 


0.098 


0.21 


n 


0.490 


0.477 


0.460 


0.456 


0.414 



TABLE I: The total (E), potential (V) and kinetic energy 
(T) per particle in Ry, and the contact value of the pair cor- 
relation function g(0), all extrapolated to the thermodynamic 
limit from unbiased RMC calculations with backflow (BF) 
nodes. We further give parameters of the momentum distri- 
bution at small k (no, and n 2 ) n(k — > 0) = no — n 2 (k / Uf) 2 , 
and at &f: n = [n(kp) + n(kp)]/2. 



Within Variational Monte Carlo (VMC), the ground 
state wave function is approximated by a trial wave 
function, "Fj^R), whereas within projector Monte Carlo 
methods, e.g. Diffusion Monte Carlo (DMC) or reptation 
Monte Carlo (RMC), the trial state is improved using 
* p oc exp[— fiH]^T] this converges exponentially fast to 
the true ground state for increasing projection time /3. 
To circumvent the so-called Fermion sign problem, cal- 
culations are done within the fixed-node approximation, 
introducing small systematic deviations from the exact 
Fermion ground state (TS]. Whenever the (approximate) 
nodes of the system are described by a determinant of 
single particle orbitals 4> n (r), the (fixed-node) ground 
state wave function, Wjv(R), of N particles at positions 
R = {r - ;}, can be written as 



* N {R) = D N exp[-U 



N\ 



Dn = det 4> n (r,+ViW}v)(l) 

nl 



where Wjv and Un are generalized backflow and Jastrow 
potentials [TB] respectively. 

From an approximate ground state wavefunction, 
'Fjv(R), we obtain the reduced single particle density ma- 
trix IP71 



Mr) = (F(R;r)> 



N ■ 
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^w(R:^-f 

*jv(R) 



(2) 



where R : r j + r indicates that the position of particle i is 



displaced by r, and 



IN 



= JdR...\^ N \ 2 /Q withQ = 



JdRI^Tvl 2 playing the role of a partition function. The 
Fourier transform of /jv(r) directly yields the momentum 



distribution, 



7 W 



of the electrons per spin 



1 

2V 



dre 



7;v(r) 



(3) 



where V is the volume. 

The large variance of the estimator of the off-diagonal 
density matrix, Eq. ([2]), makes precise calculations very 



r s 
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2 


3.99 


5 


10 


BF-RMC 


0.84(2) 


0.77(1) 


0.64(1) 


0.58(1) 


0.40(1) 


SJ-VMC 


0.894(9) 


0.82(1) 


0.69(1) 


0.61(2) 


0.45(1) 


BF-VMC 


0.86(1) 


0.78(1) 


0.65(1) 


0.59(1) 


0.41(1) 


Go Wo [25] 


0.859 


0.768 


0.646* 


0.602 


0.45 


GW [H] 




0.804 


0.702* 






gw [13 




0.846 


0.793* 






Lam [28] 


0.896 


0.814 


0.615* 


0.472 




RPA[28] 


0.843 


0.700 


0.442* 


0.323 




SJ-DMC [B] 


0.952 


0.889 




0.725 


0.593 



TABLE II: Renormalization factor, Z, extrapolated to the 
thermodynamic limit from unbiased RMC calculations with 
backflow nodes (BF-RMC), together with SJ-VMC, and BF- 
VMC results, compared with perturbative results from liter- 
ature (literature values* are at r s — 4 instead of r 3 = 3.99). 
Previous SJ-DMC results [6j used mixed estimators without 
thermodynamic limit extrapolation. 



time-consuming. To reduce the variance for homoge- 
neous systems with plane wave orbitals: </> n (r) oc e lk "' r , 
we separate the ideal gas density matrix, /i<j(r) = 
EnC( r )</ > n(0)/E„ |0„(O)| 2 , based on the estimator 



F id (R;r) = -£ 



1 v ^Z? 7V (R:r l +r;TV Ar (R)) 



D N (R;W N (R)) 



(4) 



where the determinants on the r.h.s. of Eq. Q are eval- 
uated using the backflow coordinates, Wn(R), of the di- 
agonal configuration R with un-displaced particle coor- 
dinates. Expanding it around r = 0, we can explicitly 
verify that f ld (v) = (F id (R; r)) N , so that the F — F id is a 
reduced variance estimator [T5] of the difference: Jn — fid- 
There is a problem with projecting methods to calcu- 
late properties other than the energy. Forward walking 
or reweighting methods based on using in Eq. (pj), 
become very inefficient for long projection time, since 
the variance increases exponentially with /3. To avoid 
this problem, mixed estimators, based on "F^oi are fre- 
quently used but they can introduce a systematic bias. 
Unbiased estimators for the pair correlation function, 
potential and kinetic energy have been obtained within 
RMC 12J. Based on a generalized partition function, 
Q, we extend RMC to include sampling of off-diagonal 
matrix elements [10] 



Q = 



s 



dR 1*0/3(11) | 





T,fyJ ^/^l*/MR)*r(R:r i + r)| (5) 



where s is a parameter used to optimize the efficiency 
(s = corresponds to the usual diagonal RMC [T2]). 
Similar to the worm-algorithm used in continuous Path- 
integral calculations [11] , our calculations include moves 



3 



which "open" (or "close") a path from diagonal space 
R to off-diagonal space (R, + r) . Such moves are in- 
cluded at r = and "propagated" by reptation moves 
[121 HH] to the interior of the path (r > 0). In contrast 
to previous calculations using so-called mixed estimators 
[5] , this generalization gives an unbiased estimator of the 
off-diagonal density matrix, /jv(r), and the momentum 
distribution, . Reduction of the variance based on the 
considerations above, Eq. Q, is still possible, but less 
effective. 

Quantum Monte Carlo results are obtained for typ- 
ically N < 10 3 electrons. The extrapolation to the 
thermodynamic limit introduces important quantitative 
and qualitative changes of the momentum distribution 
around the Fermi surface, Uf [5]. For a homogeneous 
periodic system, the orbitals are plane waves: (/> n (r) = 
exp[i(k„ + 9) ■ r], in the Slater determinant of Eq. QlJ> , 



where k, e G 



{( 



ni,n 2 ,n 3 



)2ttV~ 1/3 } with integer 



rii, 



and 9 can be chosen to introduce twisted boundary con- 
ditions [5) [20]. For a normal Fermi liquid, we further have 
|kj +6*| < kpi and the generalized backflow and Jastrow 
potential Wn and Un can be written exclusively in terms 
of collective coordinates pt = J2 n e r " and their gradi- 
ents [7J [TO] . Using the wavefunction "potentials" , Wn 
and Un, expressed as continuous functions in terms of 
the collective coordinates, the relation between the wave 
function in the limit N — > oo to a finite system is well de- 
fined, as it just amounts to evaluations on a denser grid 
in k-space [HI [0J ■ 

Let us first discuss the finite size scaling for a Slater- 
Jastrow (SJ) wave function: a determinant with Wn = 0, 
together with a two-body Jastrow correlation, Un — 
X)fe u kPkP-k/2V . We further assume that the func- 
tion Uk is analytically given. In our SJ-VMC calcu- 
lations, we use the Gaskell form 2nuf, J = — S' " 1 (fc) + 

1/2 

\S^ 2 {k) + 2nvk/£k] where So(k) is the ideal gas struc- 
ture factor, v k = 4ire 2 /k 2 , and e k = h 2 k 2 /2m [2Tjl22]. 
Neglecting mode-coupling between single particle modes 
in Dn and collective modes described by Un, the single 
particle density matrix, Eq. ([2]), can be approximated as 



Mr) « /c(r) 



N 



-(U' N -U N ) 



N 



(6) 



where the prime indicates the off-diagonal configuration, 
e.g. D' N = Dn(R : ri + r). Within the cumulant and 
rotating wave approximation, we then obtain an explicit 
expression, 



/cO) =s fid(r)exp{-x N (r)] (7) 
x »( r ) = V E [uk(S k -l) + nutS k ] [e ik ' r -l] (8) 



|k|<fc c 



where Sk = (picP-k) n /N is the structure factor, /jd(r) = 
2J2k<k F e ' k r /^V is the single particle density matrix of 



the corresponding ideal gas, and we have neglected con- 

1 /2 

tributions of short wave length modes, k c w 0.48r s ' kp 
PB] , Further, we can use Sk ~ [2nu k + f/S'o(fc)| _1 to 
express Sk in terms of Uk and So(k), which is based 
on assuming gaussian statistics for pk, so than Eq. ^ 
gives an explicit expression for /jv(r) ~ /c( r ) m terms 
of a given Jastrow factor. Whereas the resulting model, 
Eq. Q, depends weakly on k c , so that /jv(r) and n k 
are only qualitatively described, the size-extrapolation is 
quantitatively correct, as it is dominated by the Jastrow 
singularity u k -> (v k /2ne k ) 1/2 and S k -> {2nv k / s k )~ 1/2 
for k —¥ stemming from the plasmon contributions. 

Since we expect that mode-coupling is negligible in the 
long wave length limit, the cumulant expression, Eq. d7|), 
can be used to determine the size corrections of QMC 
calculations of the finite system 



/oo(r) 



d 3 k 



3 n k e 



-(a;oo(r) — a;jv(i")) 



(9) 



Here is the momentum distribution of the N electron 
system, defined for all values of k in a grand canonical 
ensemble using twisted boundary conditions [Hj. From 
the Fourier transform of /oo(r), Eq. we obtain the 
extrapolated momentum distribution, A related lin- 
earized expression has been used to extrapolate the mo- 
mentum distribution of the two-dimensional electron gas 
in Ref. 0. 

Following the analysis of Ref. [H] , leading order correc- 



tions to the renormalization factor, Zn 
are given by 



,N 



l k F +' 



Zoo 



Z N e-xp[-A N 
* /L d 3 q u. 



(10) 



■k/l (2tt) 3 2 

1/3 



} [1 + O {[2nu q S {q)]- 1 )] 



where c ~ 1.221 is a numerical factor to account for the 
cubic integration volume [24]. Whereas the asymptotic 
region is only reached for large systems with N 1 ^ 3 ^ 2 ^> 
1, the extrapolation based on the full expression, Eq. d9j), 
includes corrections beyond the leading order term. An- 
alyzing Eq. Q around kp, we obtain the exact leading 
order behavior with an infinite slope at kp 



n(k 



Zoo 

27 



n(k F ) 



9tt s 





' k 




log 
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i'Yi 




- 1 






kp 






kp 



(11) 



Size extrapolation, discussed above, requires the 
knowledge of the structure factor, Sk, and the Jastrow 
potential, u k , in Eq. ([8]). The QMC calculation of the 
A^-particle system allows us only to determine them on 
a finite grid in k space, but the analytic continuation to 



4 



the dense grid can be done by interpolation from their 
known behavior at small k [8]. Whereas Sk can be cal- 
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culated directly, Uk 



,SJ 



is only known explicitly for 



VMC calculations using a Slater- Jastrow trial function. 
In general, imaginary time projection and backflow in- 
troduce an effective Jastrow potential, Uk , different from 
the explicitly given form of the underlying trial wave- 
function. Expecting small changes at long wave length, 
Uk = Uf[ J + Suk, we obtain the modifications 8uk from 
from changes in the structure factor 5Sk — Sk — S§ J 
by linear response. For our purpose, mode coupling can 
be neglected, as well as deviations from gaussian statis- 
tics, so that SSk/Suk' — — ZnS^ 2 5k^> for k — > 0. There- 
fore, the effective Jastrow factor of wave functions in- 
cluding backflow and projection can be determined from 
the structure factor. 

Using SJ-VMC calculations with wf J for N = 54 
to N — 1024 electrons, we have checked that size ex- 
trapolations based on Eq. ([9| with N — 54 are reli- 
able. Thus, the more expensive backflow VMC and RMC 
calculations based on the analytical wave functions in 
Ref. [7] are only done with that size. Extrapolated re- 
sults on the total energy E, unbiased estimators from 
reptation for the potential (V) and kinetic energies (T), 
and the contact value of the pair correlation function, 
g(0), are given in table [TJ The momentum distribution 
is shown in Fig. [T] The values for the renormaliza- 
tion factor, Z, together with different perturbative re- 
sults from the literature are given in table [TTJ Table [I] 
also contains the values of the momentum distribution at 
the origin, rio, the negative slope at the origin, n^, and 



+ n k +)/2. These values can be used to param- 



eterize the momentum distribution along the lines given 
in Ref. [30], together with Z, the exact large k asymp- 
totics 29J, n(k -> oo) = (9/2)r 2 s g(0)/k 8 , and the exact 



behavior close to the Fermi surface, Eq. (111. Whereas 



the mixed estimator usually employed in DMC calcula- 
tions, introduces a small bias in the momentum distribu- 
tion, size extrapolation introduces large systematic mod- 
ifications which limit the precision of the calculations. 
Previous DMC results [5J, using mixed estimators and 
SJ nodes, suffer from these strong finite size effects and 
overestimate Z by a large amount. 

In summary, we have calcuated the momentum distri- 
bution using a new unbiased and much more accurate 
Monte Carlo method, and extrapolated the results to 
the thermodynamic limit. In particular, our data allows 
a quantitative comparison of the renormalization factor, 
Z, with approximate calculations (see table Ej). The ex- 
cellent agreement of our results with GnWn [251 |3"T1 13"2] 
over the whole metallic density region r s < 5, strongly 
indicates that vertex corrections and self-consistency is- 
sues - neither is included in GqWq — are canceling each 
other, at least close to the Fermi surface. 
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